Fourier transforms and audio

Phys481 Week 9b

import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import scipy.io.wavfile
from IPython.display import Audio
datarootdir = r'../../data/'

datafile = 'saucers-dynamics-explorer-gurnett.wav'
datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile)
S, freqs, bins, im = plt.specgram(wavedata, NFFT=1024, noverlap=512) # , Fs=fs
# Plot a spectrogram
plt.xlabel('Time') ; plt.ylabel('Frequency') ; plt.title(datafile)
#S, freqs, bins, im = plt.specgram(wavedata, NFFT=1024, Fs=fs, noverlap=512)
Audio( wavedata.T, rate= 20000)

Sinusoid with a noisy background.

Add random (“white”) noise to a sinuosid to make a signal that looks messy in the time domain. A Fourier transform will re-arrange the information, placing the harmonic energy in a few bins while the noise will be approximately uniformly distributed.

nsamples = 4096 #1024
signal = np.sin( 0.25*np.arange(0, nsamples) )
noise = np.random.randn(nsamples)

fig, axes = plt.subplots(1, 3, figsize=(15,4), sharex=True, sharey=True)

plt.sca( axes[0] )
plt.plot(signal)
plt.title('Signal')
plt.xlabel( 'Sample #' )

plt.sca( axes[1] )
plt.plot( noise )
plt.title('Noise')

plt.sca( axes[2] )
plt.plot( signal + noise )
plt.title( 'Signal\n+\nNoise' )
plt.xlim(0, 512)

(0.0, 512.0)

Fourier transform

fft_signal = np.fft.fft( signal )
fft_noise = np.fft.fft( noise )
fft_freq = np.fft.fftfreq( len(signal) )
fig, axes = plt.subplots(1, 3, figsize=(15,4), sharex=True, sharey=True)

plt.sca( axes[0] )
plt.plot( fft_freq, np.abs(fft_signal) )
plt.title('Signal')
plt.xlabel( 'frequency' )

plt.sca( axes[1] )
plt.plot( fft_freq, np.abs(fft_noise) )
plt.title('Noise')

plt.sca( axes[2] )
plt.plot( fft_freq, np.abs(fft_signal + fft_noise) )
plt.title( 'Signal\n+\nNoise' )

Text(0.5, 1.0, 'Signal\n+\nNoise')

Low-pass filter

Chop out a lot of the noise, apply the inverse Fourier transform, and get a somewhat nicer sinusoid.

fft_lowpass = np.abs(fft_freq) <= 0.075
fig, axes = plt.subplots(1, 2, figsize=(15,4))

plt.sca( axes[0] )
plt.plot( fft_freq, np.abs(fft_signal + fft_noise) )
plt.plot( fft_freq, fft_lowpass*2200, lw=3, alpha=0.5 )
plt.title( 'low-pass filter' )
plt.xlabel( 'frequency' )

plt.sca( axes[1] )
fft_lowpass_filtered = ( fft_signal + fft_noise) * fft_lowpass
lowpass_filtered = np.fft.ifft( fft_lowpass_filtered ).real
plt.plot( lowpass_filtered )
plt.xlabel( 'sample #')
plt.title( 'low-pass filtered\nSignal+Noise' )
plt.xlim(0, 1024)

(0.0, 1024.0)

import IPython.display
from ipywidgets import Layout
from IPython.display import Audio

#
#
IPython.display.display( Audio(signal, rate=4000) )
IPython.display.display( Audio(noise, rate=4000) )
IPython.display.display( Audio(signal+noise, rate=4000) )
IPython.display.display( Audio(lowpass_filtered, rate=4000) )

Audio processing

from IPython.display import Audio

# constant frequency tone
#
Audio(np.sin(np.linspace(0, 3000, 20000)), rate=20000)

# mixture of two sinusoids
#
s1 = np.sin(np.linspace(0, 3000, 20000))
s2 = np.sin(np.linspace(1000, 8000, 20000))
Audio(s1+s2, rate=20000)

# white noise - WARNING! turn down volume
#
data = np.random.uniform(-1,1,44100) # 44100 random samples between -1 and 1
scaled = np.int16(data/np.max(np.abs(data)) * 32767)
Audio( scaled, rate= 20000)

The music of the spheres?

import scipy.io.wavfile
#http://www.astrosurf.com/luxorion/audiofiles-geomagnetosphere.htm
#url = r'http://www.astrosurf.com/luxorion/Radio/chorus-cluster2-gurnett.wav'
datafile =  'chorus-cluster2-gurnett.wav'
datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile)
Audio( wavedata, rate=datarate)

# show time series
#
plt.clf() ; plt.plot( wavedata )
plt.xlabel('sample #') ; plt.ylabel('7-bit signed amplitude') ; plt.title(datafile)

Text(0.5, 1.0, 'chorus-cluster2-gurnett.wav')

# show power spectrum
#
freqdata = np.fft.fft(wavedata)
freqscale = np.fft.fftfreq( len(wavedata))
freqdata[0] = 0  # there is often a spike at zero frequence
plt.plot( freqscale, np.abs(freqdata) )
plt.xlabel('frequency') ; plt.ylabel('amplitude') ; plt.title(datafile)

print(freqdata.shape, freqscale.shape )

(2742440,) (2742440,)

# isolate low frequencies
#
lowpass = np.abs(freqscale) < 0.05
lowdata = np.real( np.fft.fft( lowpass * freqdata ) )
lowdata /= np.max(np.abs(lowdata)) / 32767
lowdata -= np.mean(lowdata)
Audio( lowdata, rate= 44100)

# isolate high frequencies
#
highpass = np.abs(freqscale) > 0.05
highdata = np.real( np.fft.fft( highpass * freqdata ) )
highdata /= np.max(np.abs(highdata)) / 32767
highdata -= np.mean(highdata)
Audio( highdata, rate= 44100)

# take a look at a small segment to see high/low frequency oscillations
plt.plot( lowdata[0:999])
plt.plot( highdata[0:999])
plt.xlabel( 'sample #' ) ; plt.title(datafile)

Text(0.5, 1.0, 'chorus-cluster2-gurnett.wav')

# 50% of low frequencies in one ear, 200% of high freqencies in the other.
#
stereodata = np.array( [0.5*lowdata, 2.0*highdata] )
stereodata.shape
Audio( stereodata, rate= 44100)

Time-frequency spectrograph

nblock = 1024
grid = np.ndarray([ len(wavedata)//nblock, nblock-1], dtype=np.float64)
for indx in range(len(grid)):
    spect = np.fft.fft( wavedata[nblock*indx : nblock*(indx+1)])
    grid[indx,:] = np.abs( spect[1:] )

plt.imshow( np.log(grid.T+0.01), aspect='auto')
plt.xlabel('time') ; plt.ylabel('frequency') ; plt.title(datafile)

Text(0.5, 1.0, 'chorus-cluster2-gurnett.wav')

# Appears to be the same plot as provided by a built-in function

fs = 44100   # sampling frequency
S, freqs, bins, im = plt.specgram(wavedata, NFFT=1024, Fs=fs, noverlap=512)

# Plot a spectrogram
plt.xlabel('Time')
plt.ylabel('Frequency')

Text(0, 0.5, 'Frequency')

datafile = 'whistler-cluster3-4-gurnett.wav'
datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile)
S, freqs, bins, im = plt.specgram(wavedata[:,0], NFFT=1024, Fs=fs, noverlap=512)

# Plot a spectrogram
plt.xlabel('Time') ; plt.ylabel('Frequency') ; plt.title(datafile)
#S, freqs, bins, im = plt.specgram(wavedata, NFFT=1024, Fs=fs, noverlap=512)

#print( wavedata.shape )
Audio( wavedata.T, rate= 20000)

# http://freewavesamples.com/alesis-fusion-fretless-bass-c3
datafile= 'Alesis-Fusion-Fretless-Bass-C3.wav'
datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile )
plt.xlabel('time') ; plt.ylabel('frequency') ; plt.title(datafile)

#wavedata.shape
S, freqs, bins, im = plt.specgram(wavedata[:,0], NFFT=1024, Fs=fs, noverlap=512)
Audio( wavedata.T, rate=datarate)

<ipython-input-52-72c780cfd596>:3: WavFileWarning: Chunk (non-data) not understood, skipping it.
  datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile )

# http://freewavesamples.com/alesis-fusion-violin-c5
datafile = 'Alesis-Fusion-Violin-C5.wav'
datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile)
wavedata.shape
S, freqs, bins, im = plt.specgram(wavedata[:,0], NFFT=1024, Fs=fs, noverlap=512)
plt.xlabel('time') ; plt.ylabel('frequency') ; plt.title(datafile)
Audio( wavedata.T, rate=datarate)

<ipython-input-53-26ce3b6b4bf7>:3: WavFileWarning: Chunk (non-data) not understood, skipping it.
  datarate, wavedata = scipy.io.wavfile.read(datarootdir + datafile)

FFT

The fast fourier transform takes a sequence of numbers and returns another sequence of the same length. Both input and output may be complex.

# FFT of 8 zeros is 8 zeros.
np.fft.fft( np.zeros(8) )

array([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

# FFT of 8 ones is a spike at zero (DC offset)
np.fft.fft( np.ones(8) )

array([8.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j])

# FFT of a DC offset is 8 ones
np.fft.fft( np.array([1,0,0,0, 0,0,0,0]) )

array([1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# FFT of a spike is a complex sinusoid
np.fft.fft( np.array([0,1,0,0, 0,0,0,0]) )

array([ 1.        +0.j        ,  0.70710678-0.70710678j,
        0.        -1.j        , -0.70710678-0.70710678j,
       -1.        +0.j        , -0.70710678+0.70710678j,
        0.        +1.j        ,  0.70710678+0.70710678j])

Formatted output

By default, the way that Numpy prints an array of complex numbers is acceptable. However, what if we wanted to make some formatting changes?

# Print one value
#
z = np.fft.fft( np.array([0,1,0,0, 0,0,0,0]) )
'{0[0]:6.2f}'.format( z )

'1.00+0.00j'

# only prints one value...
#
'{0:6.2f}'.format( *z )

'1.00+0.00j'

# only one value
#
'{:6.2f}'.format( *z )

'1.00+0.00j'

# two values
'{:6.2f} {:6.2f}'.format( *z )

'1.00+0.00j 0.71-0.71j'

# 8 values
'{:6.2f} '*8

'{:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} '

# 8 values in one string
('{:6.2f} '*8).format( *z )

'1.00+0.00j 0.71-0.71j 0.00-1.00j -0.71-0.71j -1.00+0.00j -0.71+0.71j 0.00+1.00j 0.71+0.71j '

# 8 values in 8 strings
['{:6.2f}'.format(x) for x in z]

['1.00+0.00j',
 '0.71-0.71j',
 '0.00-1.00j',
 '-0.71-0.71j',
 '-1.00+0.00j',
 '-0.71+0.71j',
 '0.00+1.00j',
 '0.71+0.71j']

' '.join( map( lambda x: '{:6.2f}'.format(x), z) )

'1.00+0.00j 0.71-0.71j 0.00-1.00j -0.71-0.71j -1.00+0.00j -0.71+0.71j 0.00+1.00j 0.71+0.71j'

Challenge: find some alternative(s)

My FORTRAN mind thinks that there must be some better way to write this. Try and find shorter/clearer methods to format a numpy array.

Back to the FFT

('{:6.2f} '*8).format( *np.fft.fft( np.array([0,1,0,0, 0,0,0,0]) ) )

'1.00+0.00j 0.71-0.71j 0.00-1.00j -0.71-0.71j -1.00+0.00j -0.71+0.71j 0.00+1.00j 0.71+0.71j '

z = np.fft.fft( np.array([0,1,0,0, 0,0,0,0]) )
plt.plot( z.real, 'go-' )
plt.plot( z.imag, 'bs--')

[<matplotlib.lines.Line2D at 0x27c8522bd08>]

z = np.fft.fft( np.array([0,0,1,0, 0,0,0,0]) )
plt.plot( z.real, 'go-' )
plt.plot( z.imag, 'bs--')

[<matplotlib.lines.Line2D at 0x27c85287948>]

z = np.fft.fft( np.array([0,0,0,1, 0,0,0,0]) )
plt.plot( z.real, 'go-' )
plt.plot( z.imag, 'bs--')

[<matplotlib.lines.Line2D at 0x27c852fad88>]

Make a set of 8-element vectors that each has 1 in a different element, and zeros in all the others. Then take the FFT of each one to get a family of corresponding frequency spectra.

x = np.array([1,0,0,0, 0,0,0,0])
for i in range(len(x)):
    z = np.roll(x,i)
    print( ('{:6.2f} '*8).format( *z ) )

  1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
  0.00   1.00   0.00   0.00   0.00   0.00   0.00   0.00
  0.00   0.00   1.00   0.00   0.00   0.00   0.00   0.00
  0.00   0.00   0.00   1.00   0.00   0.00   0.00   0.00
  0.00   0.00   0.00   0.00   1.00   0.00   0.00   0.00
  0.00   0.00   0.00   0.00   0.00   1.00   0.00   0.00
  0.00   0.00   0.00   0.00   0.00   0.00   1.00   0.00
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   1.00

x = np.array([1,0,0,0, 0,0,0,0])
for i in range(len(x)):
    z = np.roll(x,i)
    print( ('{:+6.2f}  '*8).format( *np.fft.fft(z) ) )

+1.00+0.00j  +1.00+0.00j  +1.00+0.00j  +1.00+0.00j  +1.00+0.00j  +1.00+0.00j  +1.00+0.00j  +1.00+0.00j
+1.00+0.00j  +0.71-0.71j  +0.00-1.00j  -0.71-0.71j  -1.00+0.00j  -0.71+0.71j  +0.00+1.00j  +0.71+0.71j
+1.00+0.00j  +0.00-1.00j  -1.00+0.00j  +0.00+1.00j  +1.00+0.00j  +0.00-1.00j  -1.00+0.00j  +0.00+1.00j
+1.00+0.00j  -0.71-0.71j  +0.00+1.00j  +0.71-0.71j  -1.00+0.00j  +0.71+0.71j  +0.00-1.00j  -0.71+0.71j
+1.00+0.00j  -1.00+0.00j  +1.00+0.00j  -1.00+0.00j  +1.00+0.00j  -1.00+0.00j  +1.00+0.00j  -1.00+0.00j
+1.00+0.00j  -0.71+0.71j  +0.00-1.00j  +0.71+0.71j  -1.00+0.00j  +0.71-0.71j  +0.00+1.00j  -0.71-0.71j
+1.00+0.00j  +0.00+1.00j  -1.00+0.00j  +0.00-1.00j  +1.00+0.00j  +0.00+1.00j  -1.00+0.00j  +0.00-1.00j
+1.00+0.00j  +0.71+0.71j  +0.00+1.00j  -0.71+0.71j  -1.00+0.00j  -0.71-0.71j  +0.00-1.00j  +0.71-0.71j